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Abstract 

The BFKL pomeron in perturbative QCD is plagued by the lack of unitarity and 
diffusion into the infra-red region of gluon virtualities. These two problems are inti- 
mately related. We perform numerical studies of the evolution equation proposed by 
Balitsky and Kovchegov which unitarizes the BFKL pomeron. We show how diffusion 
into the infra-red region is suppressed due to the emergence of a saturation scale and 
scaling behaviour. We study universality of this phenomenon as well as its depen- 
dence on subleading corrections to the BFKL pomeron such as the running coupling 
and kinematic constraint. These corrections are very important for phenomenological 
applications. 



1 Introduction 



The Balitsky-Fadin-Kuraev-Lipatov (BFKL) equation ^ describes the behaviour of scat- 
tering amphtudes in the hmit of high center-of-mass energy y^, in the kinematic range 
where perturbative QCD is apphcable. In the leading logarithmic approximation (LLA), 
as ^ 1 and as logs ~ 1, the behaviour is determined by the leading singularity in the 
complex angular momentum plane at j = 1 + 4A^cas In 2/7r, which corresponds to the 
vacuum quantum number exchange. In the context of QCD such exchange is called the 
BFKL pomeron. For s — > oo, the corresponding total cross section rises with energy, 
cr*°* ~ gi-i = gO-S for as ~ 0.2, violating the Froissart bound reflecting the fundamental 
principle of unitarity of the S'-matrix: 0"*°* < Clog(s). The next-to-leading corrections to 
the BFKL equation [Q] reduce the rise of the total cross section but the power-like form 
in s still holds. 

Another cumbersome property of the BFKL equation is the diffusion of gluon virtu- 
alities (dominated by the transverse momentum) into the infra-red region. The BFKL 
pomeron is made up of two interacting reggeized gluons, diagrammatically represented by 
a set of gluon ladders. With increasing energy, the distribution of the gluon virtualities 
along the ladder quickly enters the low momentum region around Aqcd with the char- 
acteristic diffusion pattern. Thus the description which is rooted in perturbative QCD is 
not consistent in the high energy limit. 

The intuitive physical picture of new phenomena, called screening or saturation effects, 
which are in agreement with the unitarity bound at high energy was provided in |p. In 
the context of the BFKL diffusion the picture looks as follows. The transverse size of 
gluons with transverse momentum k± is proportional to l/k\. For large k_\_, the BFKL 
mechanism of gluon radiation, g gg, populates the transverse space with large number 
of small size gluons per unit of rapidity Y ~ In s. The same mechanism also applies to large 
size gluons with small transverse momenta. In this case, however, the BFKL approach 
is incomplete since the gluons strongly overlap and fusion processes, gg g, are equally 
important. After taking these processes into account, unitarity is restored by taming the 
rise of the gluon density with transverse momenta below some energy-dependent scale 
Qs{Y), called saturation scale. The emergence of such scale is a fundamental property of 
the parton saturation phenomenon. The saturation scale rises monotonically with s, and 
for large enough rapidity Y for which Qs{Y) ^ ^qcd, the presented approach to the high 
energy scattering remains consistent. In this sense the solution to the unitarity problem 
cures in turn the infra-red diffusion problem. 

The parton saturation idea is realized with the help of non-linear evolution equations 
in which the gluon splitting is described by a linear term while a negative non-linear 
term results from the competing gluon recombination. The first equation of this type for 
the gluon density was postulated in |^ in the double logarithmic approximation. More 
rigorous derivation was given in . The attempts to go beyond the double-log asymptotics 
resulted in formulations based on semi-classical QCD methods ^ since the region of 
parton saturation is characterized by large values of colour fields. A different approach to 
unitary generalization of the BFKL equation was proposed by Balitsky 0, ^. By using 
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the operator product expansion for high energy scattering in QCD, he derived an infinite 
hierarchy of coupled equations for n-point Wilson-hne operators. Weigert managed to 
simphfy the form of these equations by writing them as a functional evolution equation 
for the generating functional of the Wilson-line operators |0] . The connection between the 
effective theory for the Colour Glass Condensate |^ and the evolution equation found 
by Weigert has been established in Another derivation of the small-x non-linear 

equation has been recently given in [p^j. 

The Balitsky's equations decouple in the large N^. limit. In this limit, the equation for 
the 2-point function was independently derived by Kovchegov [|l3| in the dipole picture 
The equation unitarizes the BFKL equation by including a quadratic term determined by 
the three-pomeron vertex |15]. The properties of this equation were investigated in |16|- 
|23] supporting the picture of parton saturation. The equation introduces the saturation 
scale below which the non-linear effects lead to saturation of the gluon density. Recently, 
an extended version of these equations which include cubic terms has been derived |24]. 

In this paper we concentrate on numerical studies of the Balitsky-Kovchegov evolution 
equation [13|. We illustrate the emergence of the saturation region with the saturation 
scale and scaling behaviour of the resulting cross sections. Such scaling was found in the 
DIS data at small x [^], after being inspired by the saturation model |26|. We contrast 
the properties of the non-linear BK equation with those resulting from the linear BFKL 
equation. In particular, we show the effect of unitarization on infra-red diffusion. These 
results serve as the starting point for the analysis of the Balitsky-Kovchegov equation 
with subleading corrections given by the running QCD coupling and kinematic constraint. 
Although subleading, these corrections are very important for phenomenological applica- 
tions and we carefully examine their impact onto the universality of the parton saturation 
phenomenon. 



2 The Balitsky-Kovchegov equation 

The Balitsky-Kovchegov (BK) equation was derived for the deep inelastic scattering of 
virtual photon on a large nucleus by the resummation of multiple pomeron exchanges in 
the leading logarithmic approximation (when <C 1, asln(l/x) ~ 1 and x ~ Q'^/s) and 



in the large Nc limit ||T^, 16 1 



The physical picture of this process is provided in the rest frame of the nucleus. In the 
small X limit, the virtual photon splits into a qq pair long before the interaction with the 
nucleus. Then the qq pair interacts with the nucleus. In this case the nucleus structure 
function is given by 

F2{x,Q^) = — ^ / d^rdz{\^T{z,r,Q^)\^ + \^L{z,r,Q^)\^} a{x,r) (1) 

where ^t,l are the light-cone wave functions for the transversely and longitudinally po- 
larized virtual photon, and a is the dipole cross section describing the interaction of the 
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qq dipole with the nucleus. Here, r = ri — fq is the transverse size of the dipole (ri, tq are 
the transverse coordinates of the quark and antiquark) and z is the longitudinal fraction 
of the photon's momentum carried by the quark. The sum of the squares of the photon 
wave functions (for massless quarks) is given by 

\^T? + \^L? = ^ E e}{[z' + (1 - zf] Klier) +AQ^z\l - zf Kl{er)} (2) 

where = z{l — z) Q^, Kq^i are the Bessel-Mc Donald functions and the sum is performed 
over all active quark flavours /. The dipole cross section is related to the imaginary part 
of the forward scattering amplitude of the qq dipole on the nucleus A^(r,b, y), 

a{x,r) = 2J (fhN{r,h,Y), (3) 

where b = (ri + ro)/2 is the impact parameter of the qq dipole. 

The BK equation is an evolution equation for the amplitude A^(r, b, Y). In the deriva- 
tion [13 1 it was assumed that the qq pair develops a cascade of gluons which, in the large 



Nc limit, form a system of 3-3 dipoles. Each dipole interacts independently with nucleons 
in the nucleus via two-gluon exchange. In this way the following non-linear evolution 
equation was found 



dV 

- a,. 



^;;5^rrTr'V(-- + r',b + -,F)iV(r',b + ^,n (4) 

where Og = Ncas/n, and the linear term is determined by the BFKL kernel 

(K N) (r, b, y ) = / ^ I Nir + r', b, Y) - ^ ^^'^ iV(r, b, y ) | . (5) 

Thus in the linear approximation, when each dipole scatters only once off the nucleus, the 
BFKL equation in the dipole picture is obtained. The non-linear term in (^) takes into 
account multiple scatterings and is essentially determined by the triple pomeron vertex 
flSll in the large Nc limit. Eq. (^) unitarizes the BFKL pomeron in the sense that at x — > 
and fixed, 

F2 ~ Q2 ln(l/x) . (6) 



Thus the power-like rise with energy for the BFKL pomeron is tamed |16]. 



In the infinite momentum frame, the BK equation resums fan diagrams with the BFKL 
ladders, corresponding to the pomeron, splitting into two ladders. This type of summation 
was originally proposed by Gribov, Levin and Ryskin (GLR) Q in the double logarithmic 
approximation (DLA): ag 1 and asln{l/x)ln{Q'^/AQ^j^) ~ 1. It was shown in 
that eq. (H) reduces to the GLR equation for the integrated gluon distribution. 
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Figure 1: Kinematics of the gluons in the BFKL ladder, x and xj z are the longitudinal 
momentum fractions, and k, k' and q = k — k' are transverse momenta of the gluons. 

The BK equation was derived for a large nucleus by the summation of the leading 
contributions in powers of the atomic number of the nucleus A. Other contributions, e.g. 
those corresponding to loops built out of the BFKL ladders, are suppressed by powers of 
A. However, this equation can be applied to a nucleon as a model which hints towards 
the most important contributions to the structure function at small x. The success of the 



saturation model [26| in the description of the small-x DIS data is an indication that this 
equation might also properly identify the leading contribution for a nucleon. 

A substantial complication of the analysis of eq. (^) is related to the coupling between 
the r and b variables. In our presentation we follow the approximation made in the 
earlier analyses that the dominant contribution to F2 is given by the tranverse sizes r for 
which r <^ b < R, where R is the nucleus size. In this case the coupling between r and 
b can be neglected and eq. (§) becomes independent of b. Thus, the 6-dependence can 
be suppressed in A^(r, b, y). It can then enter only through the initial condition for the 
evolution equation. In a more refined analysis, however, the exact 6-dependence should 
be taken into account. 

For an azimuthally symmetric solution, N{r,Y) = N{r,Y), eq. (^) can be rewritten 



in momentum space in a much simpler form after performing the Fourier transform [16| 

d^r , , , N(r, Y) dr 
— exp(-zkT) 2 — = / — 

where Jq is the Bessel function. In this case the following equation is obtained 



cp{k,Y) = I —exp{-ik-r)^-^ = I - J^{kr) N{r,Y), (7) 



^^^^ =as{K'®m.y) - ascl>\k,Y), (8) 

and the action of the BFKL kernel (appropriately shifted in the space of the Mellin mo- 
ments in k'^) is given by 

f^dk'^ j k'^ct>{k',Y) - k^ct>{k,Y) , k^^{k,Y) \ 
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where now k and k' are the transverse momenta of the exchanged gluons in the BFKL 
ladder, see Fig. ||. The form (P) of the BK equation serves as the starting point of our 
analysis. 



3 Saturation scale and scaling 

In this section we present some interesting properties of the solutions of the Balitsky- 
Kovchegov equation, based on a numerical analysis of eq. (|8|). In particular, we illustrate 
how the BFKL diffusion is suppressed due to non-linear effects leading to emergence of a 
saturation scale. This issue was also studied in Refs. |17|-[p^]. 

When studying diffusion in transverse momenta, it is convenient to illustrate the results 
in terms of the function (j)(k, Y) multiplied by k. The reason is very simple. In the linear 
case, when the non-linear term in eq. (^) is neglected, the resulting BFKL equation|^ has 
the asymptotic solution which is proportional to k~^. The remaining term exhibits the 
well know diffusion pattern in the variable r] = \n{k/ko). Thus for 

k <p{k, Y) = exp a,x y , _ exp - J ' ^> , 10 

where x(0) = 41n2 and x"(0) = 28(^(3), diffusion is clearly isolated. We will contrast the 
behaviour (|lO| ) with the solution of the non-linear BK equation. 



In order to solve eq. (|8|) an initial condition at 1" = has to be specified. We start 
from the condition which is concentrated at some value r/i in the logarithmic variable r] 

k^{k,Y = 0) = (5(r/-r/i). (11) 

In Fig. ^ we show the solutions for different values of Y for the linear (BFKL) and the 
non-linear (BK) equation. The BFKL solutions illustrate the features of eq. (0), i.e the 
unlimited rise with rapidity (energy) and diffusion into the domain of large and small 
values of the transverse momentum k. For the solutions of the BK equation, however, the 
small k domain is strongly suppressed and the rise with energy is significantly bounded. 
For large values of k, when the non-linearity in the BK equation is small, the BK and 
BFKL solutions coincide. 

The impact of unitarization of the BFKL pomeron on the infra-red behaviour is also 
illustrated in Fig. ^. These three dimensional plots show the function 

^{k Y) - k(t){k,Y) ^^^^ 

{Y)qy{k„,axiY),Y)' ' ' 

where kmax^^ is the value of the transverse momentum at which k (p{k, Y) has maximum^ 
Thus, k(f){k,Y) is re-normalized to unity at A; = kmax- For the linear BFKL equation, 

^ The kernel K' in eq. has a saddle point at 7 = —1/2 in the Mellin space, thus (j) ^ (fc^)"^''^. 
^ Notice that 4'(fc, Y) is no longer a solution of the non-linear eq. (fel). 
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when the asymptotic form sets in, ^(/c, Y) is given by the second exponent in (10). In this 
case kmax is independent of Y and equal to the scale introduced by the initial condition, 

kmax(Y) = h. (13) 

For the non-linear case, ^{k, Y) moves away of the infra-red region with kmaxiX) depend- 
ing on Y. This effect is also shown in Fig. |^ by making projection of constant values of 
^{k,Y) on the {rj,Y )-pla.ne. Again, for small Y, when the non-linearity in the BK equa- 
tion is negligible, the re-normalized solutions of the BFKL and the BK equations coincide. 
With increasing Y, when the non-linear effects become important, the difference between 
them in the region of small k becomes fully visible. 

Notice the parallel straight lines of constant values of the non-linear solution in a 
certain region of the plane in Fig. ^. This means that ^{k,Y) in this region is a function 
of the combination 

^ = ln(A;Ao) - Ay = In (- ^— -) A > (14) 

VA;oexp(Ay)y 

(or e^) instead of k and Y separately. 

The parameter ^ enumerates different straight lines and it is a matter of convention 
which line corresponds to ^ = 0. By a simple rescaling ^ — > C + Co; (equivalent to the 
transformation ko /coe"^") we can achieve that the straight line part of kmax{Y) corre- 
sponds to ^ = 0. In this case the re-normalized solution is a function of only one variable 
k/kmax{Y). The solution of eq. ^ has the same scaling property, 

<l,{k,Y)=<l){k/k^ax{Y)) ^ (15) 

if and only if 

(t>{kmax{Y),Y) = const , (16) 

as a function of y. As shown in Fig. ^, this condition holds true for sufficiently large 
values of Y . In this case the saturation scale can be defined 

Qs{Y) = kmaxiY) = Qo exp(Ay) , (17) 

with the exponential dependence on rapidity governed by the value of the scaling parameter 
A. The pre-factor Qq depends on the initial condition for eq. (^) while the y-dependence 
is universal. This is illustrated in Fig. ^ where the slope of the straight lines (i.e. the value 
of A) in the scaling region is the same for two different initial distributions, i.e. that given 



by eq. (11) and the Gaussian distribution 

k (j){k, y = 0) = exp \-{7] - mf/P^] • (18) 

We have verified this statement for different values of the parameters r]i and p. 

It is tempting to claim that for ^ < 0, i.e. for k < Qs{Y), scaling holds while for 
^ > it is violated. The sharp boundary, however, does not appear in reality and the 
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transition between these two regions is smooth. This is shown in Fig. ^ where the function 
k/QsiY) (p{k, Y) is plotted as a function of the scahng variable k/QgiY) for different values 
of the rapidity Y . The scaling behaviour is represented by a common line to the left of 
the point k/Qs{Y) = 1. The slow departure from scaling for k > Qs{Y) is clearly visible. 
In this sense the line k = Qs(X) is only an approximation characterizing the position of 
the transition region in the {k, y)-plane. However, the choice based on the position in k 
of the maximum of k(j){k, Y) as a function of Y is the most natural one. 

Let us finish by recalling that the scaling property of the solution of the BK equation 
results in the scaling law for the 7*p total cross section a^p ~ F2/Q'^, 

a^*p{x,Q^) = a^*p{Q'/Ql{Y)). (19) 



Such scaling law ("geometric scaling") was found in the small-x DIS data 




4 Analytical insight 

The emergence of scaling and the saturation scale may also be investigated using the 



analytic method proposed in |27]. Let us introduce the Mellin transform of the solution 
(t>{k,Y) of eq. (|): 

^{k,uj) = I dYe-'^^(l){k,Y) (20) 



where co is the Mellin variable conjugated to the rapidity Y. We adopt the following ansatz 



proposed in ref. [^| 



4>ik,u;) = 0(u;)e(^('^)-i)*, (21) 

where 7(0;) is anomalous dimension and t = ln{k'^ /k^). By taking the Mellin transform of 
eq. (P), we find 

W - a-.x(7(^))} 4>{^) e^^^")-')* = / ^ ^^"^ - -^(^'^ e(^(--')+^(-')-2)* , (22) 
where xil) is eigenvalue of the leading BFKL kernel 

X(7) = 2V'(1) - ^(7) - ^'(l - 7) , (23) 



and ip is the Euler digamma function. The right hand side of eq. ( |22|) can be evaluated in 
the saddle point approximation |27|, which results in 



{u; - a-.x(7(^))} <^(^) e^'^^^-'^' = -as ^\u;/2) e(27(-/2)-2)t | . (24) 



In the region where the non-linear term in ( ^41) can be neglected, the anomalous di- 
mension 7(0;) fulfils the equation 

u;-a-,x(7M) = 0. (25) 
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In the double logarithmic approximation, the dominant contribution to xil) comes from 
the pole 1/7, and the solution for the anomalous dimension reads 

-f{uj) = as/uj . (26) 

On the other hand, in the saturation region, where the nonlinearity cannot be neglected, 
the following condition is obtained from eq. (p^) 



7(0;) = 27(a;/2)-l, (27) 

which is satisfied by 

7(0;) = Cuj + 1 , (28) 

where C is a constant |]27[| . 

Using ( |28| ) and (^) and performing the inverse Mellin transform in w, we find the 
solution to the BK equation in the saturation region 

= I ^^e-^mu;) = J ^'^Me-^^*^''^. (29) 
We see that the solution (p(k, Y) depends only on the following combination, 



Ct + Y = 2C 



Hk/ko) + ^Y 



(30) 



which is related in a simple way to the previously defined quantity ^, eq. (0), with the 
scaling parameter A = — 1/(2C). In this way the scaling behaviour, seen in the numerical 
analysis presented in the previous section, is found. 

The solutions to eqs. (^) and (|^), as well as their derivatives with respect to 00, 
should match at some to = uJq (see for a full discussion). This allows to fix the value 
of the unknown parameter C in (^). Prom (^) one finds 

and the same relation should be valid for the solution 7(w) of eq. ( |25| ) at u; = uoc- Thus, 
by differentiation of eq. ( p5[ ) with respect to w, we find 



dxilc) X(7c) _ , . .„„^ 

= , = OLsX{lc)- (32) 

C17 1 - 7c 

In the DLA, when x{l) = 1/7) oiie obtains 

7c = - , ujc = 2as , C = — — ^ . (33) 

' 2' c s, \ ) 

Thus the scaling parameter A = — 1/(2C) = 2'as- In the case of the BK equation with the 
full BFKL kernel, we found numerically that the scaling parameter equals 



A = 2.05 as- (34) 
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This value is very close to the DLA value and the results quoted in pO| and p3l. We 



have also checked that (34) still holds in the case of collinear approximation to the BFKL 



kernel, i.e. keeping only I/7 and 1/(1 — 7) poles in eq. (p^). 

Summarizing, the non-linear BK equation has three regimes depending on the relative 
importance of the non- linear term: (I) the linear regime (when > Q'^(Y)) is dominated 
by the conventional BFKL evolution characterized by the BFKL intercept and diffusion 
in transverse momenta. In the DLA, the equation is driven by the I/7 pole of the kernel 
eigenvalue (|23|); (II) in the saturated regime (/c^ < Ql{Y)), the non- linear effects are 
strong and the gluon density is saturated. In this case, in the linear part of the BK 
equation the 1/(1 — 7) pole of xil) is important; (III) the transition between the linear 
and saturated regimes is defined by the matching (|3^) which depends on the detailed form 
of the evolution kernel. 



5 The running coupling effects 

The inclusion of the running QCD coupling (RC) into the leading order BFKL formalism 
can be treated as a partial, phenomenological way to account for a part of important 
non- leading effects. The results of the next-to-leading logarithmic calculations 1^] may be 
invoked to motivate such approach [p8| , p^ , ^]. In the context of the infra-red diffusion, 
the RC BFKL equation exposes serious problems, in particular, large sensitivity to the 
treatment of the non-perturbative region. The main result of the forthcoming analysis of 
the BK equation with the running coupling is that in this case the ambiguities related to 
the treatment of the infra-red domain are substantially reduced. 

The BFKL equation explicitly takes into account the contribution from the momenta 
down to k"^ = 0. Certainly, this is beyond the applicability of the perturbative formal- 
ism. Moreover, for the BFKL equation with the running coupling, a purely mathematical 
difficulty arises due to the Landau pole 

with 60 = 11 — 2nj/3. Therefore, a cut-off ko for gluon virtualities has to be introduced 
below which as is frozen or set to vanish. However, the BFKL amplitude becomes strongly 
sensitive to the choice of the cut-off. In the case of as{k'^) = for k'^ < fcg, the BFKL 
kernel with running coupling has a discrete, cut-off dependent spectrum. The largest 
eigenvalue Xmax, which governs the high energy behaviour of the scattering amplitude, is 
bounded by |31|: 

3.6 ,,2n 121n2 2\ 

— as{ko) < Xmax < as(Ko). (36) 

vr vr 

From this, one concludes that the equation is unstable, i.e. the pomeron intercept may 
reach an arbitrarily large value depending on the infra-red cut-off. Besides, it was found 
that the solution (f){k, Y) of the BFKL equation with the running coupling no longer 
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exhibits the diffusion pattern typical for the fixed couphng case. Namely, it takes the 
following factorized form for sufficiently large Y and k"^ > 

fe0(A;,y) ~ exp{Ay} ^ [In (k^ /kl)Y , (37) 

with A = Xmax and v = 12/(11A) — 1, so that v is positive for moderate A and becomes 
negative for A > 12/11. Thus, diffusion into infra-red is strongly enhanced in the running 
coupling case, keeping the gluon chains in the region of small virtualities. As a result, the 
high energy behaviour is governed by the non-perturbative part of the gluon configurations. 
However, it is still possible to perform collinear factorization in which the genuine non- 
perturbative contribution may be absorbed into the input for the evolution in []33| . 

The BK equation is also of perturbative origin so its validity at low fc^ may be ques- 
tioned. However, the equation is strongly constrained by unitarity requirements. In par- 
ticular, the demand for the dipole scattering amplitude to saturate to unity for large 
dipoles (the scattering matrix in the forward direction 5 = for large dipoles) results in a 
universal form of ^{k,Y) ~ InA; for small k and fixed Y. The statement about the large r 
behaviour of the dipoles is quite general. Interestingly, the Balitsky-Kovchegov equation, 
which is rooted in perturbative QCD, ensures such behaviour of the solutions at large r. 

This conclusion remains unaltered also in the scenario with the running coupling. The 
unitarity effects are strongest in the region of small values of the gluon momenta. This 
is due to the fact that the solution to the linear part of eq. (^) is proportional to fc^^ in 
the fixed coupling case, (eq. ^), and it rises faster towards low k making the non- linear 
rescattering term largest for small k. This effect tames the rapid rise of the amplitude 
governed by the running coupling taken at the cut-off scale kQ. The evolution for low k 
is essentially frozen to the known form and it is the perturbative region of k which drives 
the evolution. This was noticed in fl^ as a supersaturation phenomenon. 

In Fig. 1^ a,b we show the results of the numerical investigation of the BK equation (^) 
with the running coupling. We considered two scenarios: the cut-off regularization of the 
running coupling, 

a^(/c2) = for k^ <kl, (38) 

and freezing of the running coupling at low scales 

as{k^) ^ as{k^ + kl). (39) 

The latter option is supported by some solid theoretical arguments [U] and the QCD 
lattice data, see e.g. ref. |]35[|. The argument fe^ is given by the virtuality of the outgoing 
exchanged gluon, see Fig. and the same as(fc^) is assumed in the linear and non- 
linear terms. A different choice of the scale would correspond to a modification of the 
BFKL kernel in the NLL approximation ||2|, which does not influence significantly the 
properties of the solution. In order to investigate the sensitivity to infra-red details, two 
values of the cut-off were used: /cq =0.1 GeV^ (a) and k^ = 1 GeV^ (b). Both scenarios 
for the running coupling lead to similar results. The contours of constant values of the 
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function ^{k,Y), eq. (|T^), in the scenario ( pq ) are plotted in the y)-plane. The input 
is given by the Gaussian (|l8|) centered at ki = 10^ GeV. 

We see in Figs. § that the solution to the BFKL equation passes through a rapid 
"tunnelling" transition (the detailed discussion of this phenomenon is presented in pO[|), 
collapsing from the region of /c ~ /ci (input scale) to /c ~ feo (cut-off scale) at certain 
rapidity. As expected, for larger Y the shape in k factors out from the ^-dependence 
for the solution of the BFKL equation 4>(k,Y). In the re-normalized solution ^(k,Y), 
eq. (|l2|), the F-dependence drops out and the lines of constant values become parallel to 
the y-axis, as seen in Fig. ^. 

The behaviour of the solution of the non-linear BK equation is different. Now, we 
distinguish three regimes for ^(k,Y): 

I. For small Y, the linear BFKL regime occurs with diffusion weighted towards low k. 

11. When Y increases, the distributions broadens and the tunnelling transition occurs. 
However, this effects is different from the collapse observed in the linear BFKL case 
since the saturation scale is generated which suppresses ^{k,Y) for small k. 



111. For high Y, geometric scaling (15) appears. The saturation scale increases with Y 



and the dynamics of evolution is governed by perturbative values of k. 

The dependence on the cut-off ko is marginal for large Y in both the saturated and 
the non-saturated range of k. The statement about the independence from the infra-red 
parameter /cq also holds for the effective intercept which governs the increase of the gluon 
density with Y for the scales much larger than the saturation scale Qs{Y). 

It is interesting to estimate how the running coupling affects the Y dependence of 
the saturation scale Qs{Y). We adopt a natural approximation that the local exponent 
A(y) = dln{Qs{Y)/A) /dY takes the form X{Y) = 2a,(Q2(y)) where A = Aqcd- The 
above form is motivated by the leading logarithmic result with the fixed coupling (^) and 
its numerical verification (jsj). Thus, we have 

dHQ,{Y)/A) ^ 12 

dY boHQs{Y)/A)' ^ ' 

with the initial condition Qs{Yq) = Qq and Yq chosen in the region where scaling sets in. 
The solution takes the form 




QsiY)=Aexp{J—{Y-Yo)+LU, Y > Yq , (41) 



where Lq = ln(Qo/^)- It follows that the local exponent X(Y) decreases with increasing 
rapidity, and X{Y) ~ 1/VY for very large Y. Such dependence is indeed seen in the 
numerical analysis. 

In Fig. ^ we illustrate scaling in the running coupling case by showing the function 
{k/Qs{Y)) (p{k,Y) plotted versus k/QsiY) for different values of rapidity. QsiY) is given 
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by formula (41) with the initial condition Qs{Y = 0) = 2 GeV. The overlapping cm'ves at 
low values of the scaling variable clearly indicate that for k < Qs{Y) scaling is satisfied to 
a very good accuracy, thus justifying our ansatz ( ]4l| ) for the saturation scale. 

Summarizing, the Balitsky-Kovchegov equation avoids the self-consistency problems 
of the BFKL equation when the running coupling is included. The sensitivity to the treat- 
ment of the infra-red region is much smaller than in the linear case due to the appearance 
of the saturation scale. In contrast to the BFKL equation with the running coupling, the 
large Y asymptotics of the dipole scattering amplitude arising from the BK equation is 
governed by gluon virtualities in the perturbative domain. 



6 The kinematic constraint effects 

The BK equation (^) has been derived in the leading logarithmic approximation. It 
has been argued that the NLL corrections should be smaller that the saturation effects 
embodied in this equation. Since the corrections from the NLL BFKL kernel become 

5 /3 

important at rapidities Ia^ll ~ l/c^s , they should be parametrically smaller than the 
unitarity corrections which enter at rapidities of the order Yjj ~ l/a, ln(l/as)- However, it 
is well known that the next-to-leading corrections |^ to the BFKL kernel are numerically 
important and even make the high energy expansion unstable. Namely, the pomeron 
intercept in the NLL order, 

UNLL ^ 41n2a,(l-6.47a,), (42) 

becomes negative already for small value of a, — 0.16 ||2|. These problems have been solved 
in [28 1 where a method which stabilizes the high energy expansion by a proper resummation 



of collinear and energy-scale-dependent terms to all orders has been proposed. As a result, 
one obtains a renormalization and resummation scheme independent equation. Similar 



approaches have also been developed in 37]. In the previous section we have studied 



the NLL corrections due to the running of the coupling a,. Here, we study different 
subleading effects, namely those coming from the change of the energy scale. 

In order to verify the importance of the NLL terms we shall modify the BK equation 
in momentum space (^ by imposing so called kinematic constraint ||38| , |39| (see fig. Q), 

k'^ < e/z (43) 

on the real emission term in the linear part of the BK equation. The origin of this 
constraint lies in the fact that in the multi-Regge kinematics the exchanged gluon virtuality 
is dominated by its transverse momentum, \k''^\ ~ \k^\. After imposing the kinematical 
constraint one arrives at the following modified BK equation 



•Y 

Y) = Mk) + as I dy 







dk'^ \ k'^^{k', y) Q{Y -y- \n{k'^/k^)) - k^^{k, y) 
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where y = \n{z/x) and Y — y = \d.{1/z), see Fig. |T[ The function (f>o{k) is an initial 
condition at Y = 0. Notice that we no longer deal with a differential equation in rapidity, 
compare eq. (S). 



Condition ([4 3D leads to the following modification of the leading eigenvalue ( |23D of the 
BFKL kernel 

X(7, u;) = 2V^(1) - ^^(7) - - 7 + , (45) 

where 7 and uj are the Mellin variables conjugated to In/c^/A^ and Y, respectively. The 
indicated shift in the second '0 function accounts for the resummation to all orders of 
subleading terms resulting from the changes of the energy scale. The obtained w-dependent 
kernel ensures that the collinear limits are always correctly reproduced |28, The new 



Pomeron intercept uj^ll stays always positive and is significantly reduced as compared to 



the leading logarithmic value ujll = 41n2as ||2^, |3^. It also depends non-linearly on a 



which is a consequence of the uj dependence of the eigenvalue 

The solution of the BK equation with the kinematic constraint also exhibits the pre- 
viously discussed scaling property. Repeating the analysis of Section ^ with the BFKL 
kernel eigenvalue ( ^5[ ) instead of (|2^), we arrived at the scaling condition in the saturation 
region with the variable (^0|). Now, the value of the scaling parameter A is different. The 
matching condition (|32[) is replaced by 



9x X ( (^X 

^ ' '^c = OsXilc^^c)- (46) 



djc 1 - 7c V duj, 

Let us consider the DL approximation in which the eigenvalue (^) is expanded in 7 and 
then in u around 7 = w = 

1 7r2 

X(7,t^)- ^i^- (47) 

7 

Using this form of the BFKL eigenvalue in the matching conditions (p6|), we find 

which should be compared to the result (^) for the BK equation with the linear kernel 
without the kinematic constraint. Now, the scaling parameter A = — 1/(2C) reads 

A = , • (49) 

This approximate result suggests that the scaling parameter is reduced in comparison with 
the standard BK equation, see eq. (|3^). Moreover, the dependence of A on a, is non-linear. 

We have solved numerically the BK equation (P) with the kinematic constraint in the 
BFKL kernel and found that indeed the scaling property approximately holds. The scaling 
parameter A is found to be reduced with respect to the value (^^ for the BK equation 
without the kinematic constraint. The non-linear dependence on a, of the form (Eol) is 
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confirmed by the numerical considerations within some 15% relative deviation, towards 
smaller values of A. 

In Fig. ^ we show the solution </>(A;, Y) of the BK equation with and without the 
kinematic constraint as a function of k for different values of y. In both cases the same 
Gaussian input ( |T8|) concentrated at /ci = 1 GeV is used. The overall pattern is similar 
in both cases, with power-like fall for the large values of k and the In^k/ko) behaviour 
for small values of k. However, the transition between the two regimes at the same 
values of rapidity occurs for smaller values of k in the case of solution with the kinematic 
constraint. This means that the saturation scale generated by the BK equation is smaller 
when the kinematic constraint is present. This is a phenomenon analogous to the well 
known reduction of the leading order BFKL pomeron intercept by including the next-to- 
leading corrections. 

Prom Fig. ^ it is seen that the kinematic constraint is potentially very important for 
the phenomenological studies based on the BK equation. Its effect manifests for instance 
in a slower increase of the saturation scale than in the standard case. The impact of 
the kinematic constraint is even more important in the linear regime, where the solution 
(j){k, Y) drops by two orders of magnitude at y = 20 (for Og = 0.2) when the constraint 
is included. However, the general picture of the saturation process remains unaltered. 
So, future reliable phenomenological applications of the BK equation to should involve 
the kinematic constraint or another representation of resummed non-leading corrections 
to the BFKL part of the evolution kernel. It might also be valuable to combine the 
kinematic constraint and the running coupling effects in the BK equation. 

7 Summary 

We presented detailed numerical analysis of the non-linear Balitsky-Kovchegov (BK) equa- 
tion. We showed how diffusion into the infra-red domain for the linear BFKL equation 
is suppressed for the BK equation due to the emergence of the saturation scale Qs(Y), 
which depends on the rapidity Y. The saturation scale divides the {k, y)-plane into a 
region governed by the linear BFKL evolution when k > Qs{Y), and the region in which 
gluon saturation occurs when k < Qs{Y). In the latter region, the solution is a function 
of only one variable k/Qs{Y), instead of k and Y separately (scaling behaviour). Our 
results are fully compatible with the previous analyses performed at the leading order. 
The scaling parameter A ~ 2as which governs the Y dependence of the saturation scale, 
Qs ~ exp (Ay), is universal and does not depend on the form of the initial conditions for 
the evolution. 

We also investigated the influence of the NLL corrections to the BFKL kernel such 
as the running coupling and kinematic constraint. Similarly to the fixed coupling case, 
when as runs, the diffusion of gluon virtualities into the infra-red domain is avoided by 
generation of the saturation scale Qs- We found, however, a different dependence of Qs on 
the rapidity, see eq. (|4l]). Thus, the non- linear evolution at high rapidity is driven by gluon 
momenta in the perturbative range independently of the initial condition for the evolution. 
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The details of the BK solutions in the running coupling case are quite insensitive to the 
regularization of Ug at small scales. 

We also studied the impact of the kinematic constraint onto the non-linear BK equa- 
tion and found that an approximate scaling property still holds. However, the kinematic 
constraint slows down the rise of the gluon density with increasing rapidity (energy). 
Thus, the saturation scale increases slower with rapidity than in the unconstrained case, 
see eq. (p^). 

We derived and tested numerically the analytic formulae describing the dependence of 
the saturation scale on rapidity when either the kinematic constraint or the running QCD 
coupling are included in the BK equation. This is the first systematic study of impact of 
those subleading corrections on the gluon density saturation in the Balitsky-Kovchegov 
framework. These results suggest that the subleading effects are very important for the 
construction of a reliable phenomenology. 
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Figure 2: The functions k(j){k, Y) constructed from solutions to the BFKL and the Balitsky- 
Kovchegov equations with the input for different values of the evolution parameter 
Y = lii(l/x) ranging from 1 to 10. The coupling constant Os = 0.2. 
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Figure 3: The re-normalized solutions "^(kjY) from the input ([TJj for the BFKL and 
the Balitsky-Kovchegov equations as a function of the rapidity Y and i] = In^k/ko) with 
ko = 10-10 GeV 
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Figure 4: The lines of constant values of the BFKL and the BK re-normalized solutions 
*(A;,y) (Y = ln{l/x)) in the {logiQ{k),logiQ{l/ x))-plane. 
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Figure 5: Scaling condition (j7^ as a function ofY. 



21 




-2 -1 012345678 



logi„(k/lGeV) 



Figure 6: The lines of constant values of the Balitsky-Kovchegov re-normalized solutions 
^{k,Y) in the {logiQ{k),logiQ{l /x))-plane for the delta-like input (11) and the Gaussian 
input i^T^). 
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Figure 7: The function (k/Qs{Y)) (/){k,Y) plotted versus k/Qs{Y) for different values of 
rapidity Y ranging from 10 to 23. The saturation scale Qs{Y) corresponds to the position 
of the maximum of the function k (p{k, Y). 
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Figure 8: The re-normalized solutions "^{k^Y) from the input ^Tl\) for the Balitsky- 
Kovchegov equation with the running coupling constant and the infra-red cut-off = 
O.lGeV^ (a) and k^ = 1 GeV^ (b) as functions of logi^il/x) and logio(fc/l GeV). 
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Figure 9: The function {k/Qs{Y)) (f){k,Y) in the running coupling case, plotted versus 
k/Qs(Y) for different values of rapidity Y ranging from 15 to 32. The saturation scale 
Qs{Y) is taken from equation with the initial condition QsiY = 0) = 2.0 GeV . 
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Figure 10: Comparison of solutions of the Balitsky-Kovchegov equation with and without 
the kinematic constraint, plotted versus k for different values of rapidity Yi = 20 + i * 2, 
where i goes from to 10. The dotted line - the input distribution atY = 0; the solid line 
- the result of the calculation with the kinematic constraint, the dashed line - the result of 
the calculation without the kinematic constraint. The coupling constant as = 0.2 is fixed. 
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